/*!
 ***************************************************************************
 * \file   convection1.cpp
 * Forced Convection-Diffusion in 1D.
 *
 ***************************************************************************
 * \author    Luis M. de la Cruz Salas [ Sat Apr 19 18:47:23 BST 2008 ]
 ***************************************************************************
 * This code solves the next equation:
 * \f[ 
 * \frac{\partial \psi}{\partial t} + 
 * u \frac{\partial \psi}{\partial x} = \frac{\partial^2 \psi}{\partial x^2}
 * \f]. 
 * \par
 * This equation is solved in the domain \f$ x \in [0,\infty] \f$. 
 * The velocity \f$ u \f$ is equal to 1.
 *
 * The equation can be solved using Upwind, Central Average (CDS) or QUICK 
 * schemes for the approximation of convective terms. CDS and QUICK can
 * be used with \c CDS_CoDi and \c QUICK_CoDi, or with \c CDS_Hay and 
 * \c QUICK_Hay. The latter two are the deferred implementation from
 * Hayase et al. 1992, which are more stable. 
 *  \par Compiling and running
 *  Modify the variables BASE and BLITZ in the file \c tuna-cfd-rules.in 
 *  according to your installation and then type the next commands:
 * \par
   \verbatim
   % make
   % ./convection1 \endverbatim
 * The results are stored in \c Data1 .
 ***************************************************************************
 */

#include <iostream>
#include <cmath>
#include "Meshes/Uniform.hpp"
#include "Equations/ScalarEquation.hpp"

#ifdef WITH_CDS_CoDi
#include "Schemes/CDS_CoDi.hpp"
#endif

#ifdef WITH_CDS_Hay
#include "Schemes/CDS_Hay.hpp"
#endif

#ifdef WITH_Upwind_CoDi
#include "Schemes/Upwind_CoDi.hpp"
#endif

#ifdef WITH_Quick_CoDi
#include "Schemes/Quick_CoDi.hpp"
#endif

#ifdef WITH_Quick_Hay
#include "Schemes/Quick_Hay.hpp"
#endif

#include "Solvers/TDMA.hpp"
#include "Utils/inout.hpp" 
#include "Utils/num_utils.hpp"
#include "Utils/GNUplot.hpp"
using namespace Tuna;
  
// Scalar field in one dimension
typedef TunaArray<double,1>::huge ScalarField1D;

void exactSolution(ScalarField1D &, double, double, double, double);

int pausa;

int main() 
{

#ifdef WITH_GNUPLOT
    GNUplot plotter;
    int stop;
    plotter("set grid");
#endif

    double length_x;
    int num_nodes_x;
    double right_wall = 0.0, left_wall = 1.0;
    double dt, error = 1.0, residual, G;
    int nmax, frequency, iteration = 1;

// Reading data from the input file
    std::ifstream input_cin ("input1"); 
    input_cin >> length_x
	      >> num_nodes_x 
	      >> G
	      >> dt
	      >> nmax
	      >> frequency;
    input_cin.close();

    StructuredMesh<Uniform<double, 1> > mesh(length_x, num_nodes_x);
    double dx = mesh.getDelta(X);
    int num_vols_x = mesh.getExtentVolumes(X);
    mesh.print();

    ScalarField1D T(mesh.getExtentVolumes());
    ScalarField1D u(mesh.getExtentVolumes());
    ScalarField1D us(num_nodes_x);         // staggered u-velocity 

    ScalarField1D Tn(mesh.getExtentNodes());
    ScalarField1D Ta(mesh.getExtentNodes());
    ScalarField1D un(mesh.getExtentNodes());

    T = 0;
    T(T.lbound(firstDim)) = left_wall;   
    T(T.ubound(firstDim)) = right_wall;

    NumUtils::interpolateToNodes(Tn, T);
    InOut::writeToFile(Tn, iteration, "./Data1/temp.", dx);

    un = 1;
    NumUtils::staggerFromNodesX(us, un);

    SparseMatrix< Diagonal<double, 1 > > A(num_nodes_x);

    ScalarField1D              b(num_nodes_x);
        
#ifdef WITH_CDS_CoDi
    ScalarEquation<CDS_CoDi<double,1> > conv_diff(T, A, b, mesh.getDeltas()); 
#endif

#ifdef WITH_CDS_Hay
    ScalarEquation<CDS_Hay<double,1> > conv_diff(T, A, b, mesh.getDeltas()); 
#endif

#ifdef WITH_Upwind_CoDi
    ScalarEquation<Upwind_CoDi<double,1> > conv_diff(T, A, b, mesh.getDeltas());
#endif

#ifdef WITH_Quick_CoDi
    ScalarEquation<Quick_CoDi<double,1> > conv_diff(T, A, b, mesh.getDeltas());    
#endif

#ifdef WITH_Quick_Hay
    ScalarEquation<Quick_Hay<double,1> > conv_diff(T, A, b, mesh.getDeltas());    
#endif

    conv_diff.setGamma(G);
    conv_diff.setDeltaTime(dt);
    conv_diff.setDirichlet(LEFT_WALL);
    conv_diff.setDirichlet(RIGHT_WALL);
    conv_diff.setUvelocity(us);
    conv_diff.print();

    ofstream error_fp("./Data1/error"), residual_fp("./Data1/residual");
    int Nerror;
    double errorL1, errorL2, errorLmax, errorL2_I, errorMS, errorMS_I;
    double AeL1, AeL2, AeLmax, ReL1, ReL2, ReLmax;
    double normL1, normL2, normLmax;
    timer time;
    double total_time = 0.0;

    while ( iteration <= nmax ) {
        time.tic();
	conv_diff.calcCoefficients();
	Solver::TDMA1D(conv_diff);
	residual = conv_diff.calcResidual();
	conv_diff.update();
	total_time += time.toc();

	if (!(iteration % frequency)) {
	  exactSolution(Ta, 1, iteration * dt , G, dx);
	  NumUtils::interpolateToNodes(Tn, T);
	  errorL1   = NumUtils::calcErrorL1(Tn, Ta);
	  errorL2   = NumUtils::calcErrorL2(Tn, Ta);
	  errorLmax = NumUtils::calcErrorLmax(Tn, Ta);
	  AeL1      = NumUtils::calcAbsoluteError_L1(Tn, Ta);
	  AeL2      = NumUtils::calcAbsoluteError_L2(Tn, Ta);
	  AeLmax    = NumUtils::calcAbsoluteError_Lmax(Tn, Ta);
	  normL1    = NumUtils::normL1(Ta);
	  normL2    = NumUtils::normL2(Ta);
	  normLmax  = NumUtils::normLmax(Ta);
	  ReL1      = AeL1 / normL1;
	  ReL2      = AeL2 / normL2;
	  ReLmax    = AeLmax / normLmax;
	  errorMS   = errorL2 * errorL2;
	  
	  /* */
	  errorL2_I = 0.0; Nerror = 0;	
	  for(int i = Tn.lbound(firstDim)+1; i <= Tn.ubound(firstDim)-1; ++i, ++Nerror) {
	    errorL2_I += (Tn(i) - Ta(i)) * (Tn(i) - Ta(i));
	  }
	  errorL2_I = sqrt(errorL2_I / Nerror);
	  errorMS_I = errorL2_I * errorL2_I;
	  /* */
	  //	error_fp << dt * iteration << "\t" << error << std::endl;
	  
	  cout << "\n ---- Iteration = " << iteration
	       << "\n\n Residual = " << residual << "\n"
	       << " Error L1 = " << errorL1 << "\n" 
	       << " Error L2 = " << errorL2 << " I = " << errorL2_I << "\n" 
	       << " Error Lmax = " << errorLmax << "\n" 
	       << " Error MS = " << errorMS << " I = " << errorMS_I << "\n"
 	       << " Error Absolute L1 = " << AeL1 << "\n" 
 	       << " Error Absolute L2 = " << AeL2 << "\n" 
 	       << " Error Absolute Lmax = " << AeLmax << "\n" 
 	       << " Error Relative L1 = " << ReL1 << "\n" 
 	       << " Error Relative L2 = " << ReL2 << "\n" 
 	       << " Error Relative Lmax = " << ReLmax << "\n"
	       << " Total time = " << total_time << "\n";
	  
	  cout << errorL1 << " " 
	       << errorL2 << " " 
	       << ReL1 << " "
	       << ReL2 << " "
	       << ReLmax << " "
	       << ReLmax << " "
	       << residual << " "
 	       << total_time << "\n\n";

	  InOut::writeToFile(Ta, 1, "./Data1/exac.", dx);
	  InOut::writeToFile(Tn, 1, "./Data1/temp.", dx);
	  
	  /* */
#ifdef WITH_GNUPLOT
	  std::cout << "\n >---> Press <1> and then <enter> to continue \n";
	  std::cin >> stop;	
	  plotter("p [0:2.5][-0.25:1.5] \"Data1/temp.00001\" w l title \"Numerical Simulation\", \"Data1/exac.00001\" title \"Exact solution\" w l ");
#endif
	  /* */
	}	    

	cout << "\n ---- Iteration = " << iteration;
	++iteration;
    }

    return 0;
}

void exactSolution(ScalarField1D &Ta, double vel, double t, double D, double dx)
{
    double x;
    for(int i = Ta.lbound(firstDim); i <= Ta.ubound(firstDim); ++i) {
      x = i * dx;
      Ta(i) = 0.5 * ( erfc( 0.5 * (x - vel * t) / sqrt(D * t) ) + 
		      exp(vel * x) * exp(-D) * 
		      erfc( 0.5 * (x + vel * t) / sqrt(D * t)) );
    }
}





